Slow dynamics under gravity: a nonlinear diffusion model 
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We present an analytical and numerical study of a nonlinear diffusion model which describes 
density relaxation of loosely packed particles under gravity and weak random (thermal) vibration, 
and compare the results with Monte Carlo simulations of a lattice gas under gravity. The dynamical 
equation can be thought of as a local density functional theory for a class of lattice gases used to 
model slow relaxation of glassy and granular materials. The theory predicts a jamming transition 
line between a low density fluid phase and a high density glassy regime, characterized by diverging 
relaxation time and logarithmic or power-law compaction according to the specific form of the 
diffusion coefficient. In particular, we show that the model exhibits history dependent properties, 
such as quasi reversible-irreversible cycle and memory effects - as observed in recent experiments, 
and dynamical heterogeneities. 



I. INTRODUCTION 

The dynamics of granular matter has received consid- 
erable attention in the past few years as it poses inter- 
esting problems from a theoretical point of view, besides 
its relevance to industrial applications Q, 0, 0, H- At 
high density, excluded-volume interactions play a crucial 
role in the formation of disordered, amorphous granular 
packings. In fact, the analogy between slowly compact- 
ing granular materials and other disordered systems like 
glasses has been early recognized 0, 0, and has moti- 
vated several experiments in which slow compaction and 
history dependence have been investigated in great de- 
tail IHIEIaliEini 

Granular and glassy systems share the important fea- 
ture of having an exponentially large number (in the sys- 
tem size) of different mechanically stable packings. Mi- 
croscopically, this property can be thought of as gener- 
ated by geometric frustration or kinetic constraints on 
the possible moves or positions of particles. This in turn 
leads, at high packing density, to a vanishing particle mo- 
bility which is the distinctive macroscopic manifestation 
of slow relaxation and jamming transition (dynamical ar- 
rest). 

In general, two ingredients are responsible for the un- 
usual behavior of a granular material. First, collisions 
between the particles are inelastic, and energy has to 
be constantly pumped into the system. Second, at high 
packing density, the excluded volume, and the associ- 
ated cage effect, is very similar to the one observed in 
structural glasses A number of schematic lattice- 
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gas models Q E E 03, 13 G3 has shown indeed 
that the main features of irreversible compaction do not 
depend, to a certain extent, on the dissipation mecha- 
nism (which is often assumed to be, for simplicity, of 
thermal nature), but can be understood solely in terms 
of steric hindrance. In this paper, we shall be concerned 
precisely with this quasi-static flow regime, which in the 
case of dense systems is the relevant one. Even with this 
simplifying assumption the detailed correspondence with 
mean-field, mode-coupling approaches of glassy dynam- 
ics remains however problematic because the presence of 
gravity leads to a non-trivial dependence on the spatial 
variable for the basic observables. Other complications 
may further arise from the presence of boundary condi- 
tions. At the present stage, coarse-grained approaches 
based on real-space diffusion equations can therefore be 
very useful in the theoretical interpretation of experimen- 
tal results and to disentagle the glassy features which 
are inherent to the compaction dynamics from the ones 
which depend on the specific energy injection/dissipation 
mechanism. 

Some typical questions in slow granular dynamics that 
one is concerned with are the: 

• origin of the logarithmic compaction law; 

• scaling behavior of the aging dynamics; 

• reversible-irreversible cycle; 

• memory phenomena; 

• and dynamical heterogeneities. 

We have recently addressed some of these issues by 
studying the dynamics of a kinetic free-volume model 
for granular media [3, and proposing an analytical ap- 
proach based on a dynamical local density functional the- 
ory |2(J. We have precisely characterized the way com- 
paction and aging depend on the particle mobility of a 
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homogeneous system. Specifically, there are two relax- 
ation regimes, fast and slow, separated by a dynamic 
jamming transition. The compaction law in the slow dy- 
namic regime depends upon the particle mobility: we 
find that logarithmic compaction and simple aging are 
intimately related to the Vogel-Fulcher law, while power- 
law compaction and superaging behavior occur in the 
presence of a power-law vanishing mobility. The objec- 
tive of this paper is to extend our previous work [jja, |2(| , 
and in particular to address other history dependent phe- 
nomena, such as reversible-irreversible cycles @, |tj lTo| , 
memory phenomena pll I22T I23L |24| and dynamical het- 
erogeneities [25l |2(| , which have recently attracted some 
attention. 

In the next section we shall introduce the nonlinear dif- 
fusion model and the related lattice-gas which is thought 
to be its microscopic realization. The stationary den- 
sity profiles, along with the jamming transition will be 
presented in section ITTT1 In section ITvl the possible sce- 
narios for the time evolution of packing density are de- 
scribed. Reversible-irreversible cycles and memory phe- 
nomena will be discussed in section while the scaling 
behavior of the aging dynamics is addressed in section lVll 
Section IVIII discusses the dynamical heterogeneities in 
presence of gravity and, finally, in section IVIIII the con- 
clusions will be presented. 



II. MODELS 
A. Nonlinear diffusion equation 

We assume that the dynamical evolution of the local 
particle density p{z, t) is governed by the continuity equa- 
tion, 



dp(z,t) dJ(z,t) 



= 



dt dz 

with the particle current J(z, t) given by the Fick's law, 

dp{z, t) 



J(z,t) = -T(p)- 



dz 



where T(p) is the Onsager mobility and p(z,t) = ^£ is 
the local chemical potential. The only interaction be- 
tween the grains we consider is the hard core repulsion, 
for which the exact lattice Hclmholtz free energy func- 
tional is I23 



PF\p(z,t)] = [ dz[ 1Z p-S{ P )] 
Jo 



(1) 



where the entropy S(p) is given by 

S(p)=-pmp-(l-p)ln(l-p). (2) 

For highly packed hard-sphere systems, theoretical and 
experimental studies suggest that the diffusion coefficient 



vanishes as a power-law |2g, |2!j, |3fJ . Hence we will as- 
sume that the mobility T(p) vanishes as 



Pc 



T(p)=T p 1- 



and remains zero for p > p . Below we will also discuss 
another possible functional form for the mobility which 
is commonly encountered in systems of particles with 
anisotropic shape (e.g. rods). Note that the above func- 
tional form of mobility has an implicit dependence on the 
height, because the density profile is typically inhomoge- 
neous (p — p{z)) due to the driving force and boundary 
conditions. The use of a local density approximation for 
the mobility will be justified by the comparison of the- 
oretical predictions with Monte Carlo simulations of a 
lattice gas model which exhibits a vanishing diffusion co- 
efficient at a threshold density p c |31| . Substituting T(p) 
into the continuity equation we are led to 




1 dp 

1 — p dz 



IP 



(3) 



where the time is now measured in units of l/T k B T. 
This equation has to be completed by specifying the 
boundary conditions. We will discuss two simple cases 
corresponding to open and closed systems. In both situ- 
ations one boundary condition requires the vanishing of 
the current at the bottom layer z = 0, J(0,t) — for 
any time t. If the top layer z — H is in contact with a 
particle reservoir at density p R (open system), the other 
boundary condition reads p{H, t) = p n for all t; while for 
a closed system in which the total number of particles is 
kept constant, the second boundary condition leads to a 
vanishing current also at z = H, J(H, t) = 0. Although 
no closed analytical solutions of Eq. © is found, it is pos- 
sible to characterize its asymptotic long time regime by 
an explicit calculation of density relaxation and two-time 
mean-square displacement. 



B. Microscopic lattice gas 

The simplest microscopic realization of the nonlinear 
diffusion equation one can imagine is provided by a lattice 
gas having a vanishing diffusion coefficient above a cer- 
tain threshold density (see j^] and references therein) . A 
paradigmatic example is the kinetically constrained lat- 
tice gas model devised by Kob and Andersen 0]. The 
model was originally introduced with the purpose to test 
the predictions of the mode-coupling theory for super- 
cooled liquids [2^] . The system consists of N particles 
on a lattice with at most one particle per site and no 
other static interactions between the particles, that is the 
Hamiltonian is Ji = 0. The microscopic dynamic is as 
follows: at each time step a particle and one of its neigh- 
boring sites are chosen at random; the particle can move 
to the new site if this site is empty and if the particle has 



3 



less than v nearest neighbors occupied before and after 
the move. This kinetic rule is time-reversible and the de- 
tailed balance is satisfied. At high densities, the dynam- 
ics slows down because the reduced free- volume makes it 
harder for a particle to satisfy the dynamic constraints. 
There exists a critical density p c above which the par- 
ticles are so interlocked that no macroscopic structural 
rearrangement is possible and the mobility falls to zero 
as a power law, D(p) ~ (p c — p)^. For the simple cubic 
lattice |3l|. p r ~ 0.88, while for the body centered cubic 
(BCC) lattice one gets p c ~ 0.84, see Fig[T] The critical 
density p c is therefore non universal, depending both on 
the lattice structure and the particular choice of the dy- 
namical constraint parameter v. In both cases, however, 
the value of the exponent is consistent with <p ~ 3.1. The 
universality of the exponent 4> was recently suggested by 
Imparato and Peliti, who studied the diffusion on a face 
centered cubic lattice for several choices of v [3^ . 

Since much of dynamical properties of both struc- 
tural glasses and dense granular materials are dictated 
by steric constraints, we have generalized the Kob- 
Andersen model in a gravitational field. The Hamilto- 
nian now is 



global density or the chemical potential, according to the 
statistical ensemble. 
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FIG. 1: The diffusion coefficient of a BCC lattice (L = 32) as 
a function of the density. The full line is a power-law fit with 
p c ~ 0.84 and <j> = 3.1. 



III. JAMMING TRANSITION AND DENSITY 
PROFILES 



where n.; = 0,1 is the occupation variable of the i-th 
site whose height is Zi, 7 = mg/k B T is the inverse grav- 
itational length and g is the constant gravitational field 
acting in the — z direction. Since we are interested in 
the slow compaction regime, the energy dissipation due 
to inelastic collisions is ignored. The thermal energy of 
the grains is negligible and T is neither the physical tem- 
perature nor the "granular temperature" usually associ- 
ated with the average kinetic energy, but rather a func- 
tion of the externally imposed vibration intensity. In 
other words, we assume that the random diffusive mo- 
tion of "grains" produced by the mechanical vibrations 
of the box can be modeled as a thermal bath of temper- 
ature T . The particles satisfying the kinetic con- 
straints can move according to the Metropolis rule with 
probability min[l,a; _AA ], where Ah = ±1 is the verti- 
cal displacement for the attempted elementary move and 
x = exp(— 7) represents the "vibration amplitude". Par- 
ticles are confined in a box closed at the bottom and with 
periodic boundary condition in the horizontal direction. 
At the top, the box can be either closed or in contact 
with a particle reservoir. We set the constraint threshold 
at v = 5. The Markov process generated by the kinetic 
rules is irreducible on the full configuration space ^| , the 
static properties of the model are those of a lattice-gas 
of non-interacting particles in a gravitational field, and 
these can be easily computed. For example, the mean 
occupation of each level is: 



p{z) = 



1 



1 + e'^+'i) 



(5) 



where the Lagrange multiplier r\ is determined by the 



A. Open system 

The stationary state of Eq. is obtained when 
dp/dt — 0, which implies that J(z, 00) = for all z. 
Imposing the stationarity condition, depending on the 
value of 7, two very distinct types of stationary profiles 
are found. For high values of 7 the system is in equi- 
librium and the profile is given by Eq. [5] As the vibra- 
tion is lowered, a homogeneous region of constant density 
p[z, 00) = p c develops below the height zq. The station- 
ary solution of Eq. |3| becomes, 



Poo(z) = < 



Pc Z < Z 
1 



1 + exp (jz + 77) ' 



Z > Zq 



(6) 



The values of zq and rj are obtained from the boundary 
conditions. If the top of the box is connected to a particle 
reservoir then Poo(H) — p R . This and the continuity of 
the density profile at z — z lead to, 



n 



H 



In- 



Iln^i 

7 Pc(l- 

- Pr 



Pc) 



Pr) 



Pr 



(7) 
(8) 



The jamming transition corresponds to the locus in the 
parameter space (7, p R ), at which the z = layer attains 
the critical density so that p(0,oo) = p c . This happens 
when 



, s 1 , Pc(l ~Pr) 
7c Pr = T7 m —71 C 

H p R (l - p c ) 



(9) 
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When 7 > 7 c (Pr) the density at the bottom of the 
box is close to the critical and dynamics becomes slug- 
gish. On the other hand, above the critical temperature, 
7 < 7c(/3r), all the layers have densities smaller than p c 
and the system easily attains equilibrium. The critical 
line Eq. ED is plotted, as a function of p R , in Fig. El No- 
tice that for the undriven case (7 = 0) the transition 
only occurs if p R = p c . It is important to stress that for 
7 > 7c (Pr) the stationary profiles Eq. 0are not equiva- 
lent to the equilibrium ones since they do not minimize 
the Helmholtz free energy functional Eq. ^ This is so 
because the system is not able to achieve, dynamically, 
densities higher than p c . In the zero gravity case 7 = 0, 
the stationary profile is flat, poo(z) = p R . 
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FIG. 2: Jamming transition line y c (multiplied by H) as a 
function of p R , from Eq.|5] This line separates regions of slow 
(glassy) and fast (fluid) relaxation. 

In Fig. we compare the stationary profiles with the 
ones found in Monte Carlo simulation of the gravity- 
driven KA model on the BCC lattice 0, |2(j at large 
times. A very good agreement is obtained with no ad- 
justable parameters. As discussed in the next section, 
the simulations were carried on connecting the topmost 
layer of the system to a reservoir of particles. 



B. Closed system 
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FIG. 3: Stationary profiles above and below the critical line 
for an open system in contact with a particle reservoir at 
density p R = 0.1. The symbols are the densities obtained from 
the simulation on a BCC lattice and the solid lines are the 
theoretical results. For 7 = 0.072 > y c , squares, the longest 
time shown is 10 6 MCS while for 7 = 0.041 < 7 C , circles, the 
time is 10 5 MCS. Notice that the simulation profile for 7 > j c 
is not yet stationary. 



This line is depicted in fig. 0] as a function of the aver- 
age density p. An example of a stationary profile for a 
fixed number of particles is shown in fig. [S] Simulations 
on bidimensional hard spheres also present a profile com- 
patible with an almost flat part plus an interface |34j |. 
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For a closed system of volume V (the height H times 
the basis area) and a fixed number of particles N = ~pV 
one proceeds in a way analogous to the previous section. 
We find that Zq and rj satisfy the coupled equations 



7 Pc 



1 

1 



rj = In ■ 



-y[pH+z (l-p c )] 



X _ e -r(PcZo-Hp) 



(10) 

(11) 



while the locus of the jamming transitions satisfies the 
implicit equation 



1 



7c 



ln(l 

pH [ 



Pc + p c e 



(12) 



FIG. 4: Jamming transition line, 7 C (multiplied by H), as 
a function of the average density p. The divergence near the 
origin scales as p _1 , different from the logarithmic divergence 
found for the system in contact with a reservoir, Eq.El Notice 
that in this case the density parameter is the average total 
density, p, while in the open case it is the reservoir density 

PR. 

We now turn to the behavior of the density profile 
found in the Monte Carlo simulation. We first let the 
system evolve at x — until a mechanically stable config- 
uration in which no particle can move down, is achieved. 
This state is metastable because at x = 0, there is a sin- 
gle ground state with equilibrium bulk packing density 
p = 1. Experimentally, it has been verified that the ini- 
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FIG. 5: Examples of density profiles above and below the 
critical line for a system closed after its density achieved p = 
0.5. The tapping amplitude in the glassy phase is 7 = 0.4, 
while in the fluid phase it is 7 = 0.03; system height H — 80. 
Notice that for 7 > 7 C the profile at finite times (10 6 MCS in 
this case) is neither stationary nor flat. 



tial density depends on the preparation procedure |35j| . 
in particular, on the box filling rate, that is, the num- 
ber of particles which fall per unit of time. If particles 
are poured one at a time, that is, they are individually 
handled, then the falling trajectories are independent. 
By doing so in the gravity driven KA model, the kinetic 
constraints would always be satisfied since the sites above 
the falling particle would be empty, and the system would 
achieve a fully compacted state. By pouring more than 
one particle at a time, their trajectories may interfere, 
preventing the system from achieving the highest possi- 
ble density. This collective handling of particles ma y b e 
implemented in several ways. For example, in ref. [15|, 
in order to avoid this highly compacted initial state, all 
the particles were poured together, by placing them in 
the upper half of the box and letting them fall randomly 
under the action of the gravitational field, until a state 
where no particles can fall further is achieved. In this 
case the initial average packing density is p r i p ~ 0.707, 
roughly corresponding to a random loose packed state. 
Once the system is prepared, the vibration at a fixed 
amplitude x is turned on. Another possibility, pointed 
out in the previous section is to use a reservoir that may 
be left open forever or be closed after a predetermined 
number of particles has entered the box. This is also 
a collective way of handling the particles, with the ad- 
vantage that the initial flux of particles can be tuned, 
being intermediate between all particles falling at once 
and one at a time. Indeed, as the system is less con- 
strained than when all the particles are falling at once, 
it is able to achieve a larger bulk density and the profile 
is closer to the stationary one (see figs. Eland EJ). More- 
over, for the same reason, the structured region at the 
bottom, discussed in |l5j , is either absent or significantly 
less pronounced. 



Near the top of the granular pile, a dense interfacial 
layer forms as a result of an increased mobility due to 
the low density of particles in the top most layers. This 
dense layer becomes more compact with time hindering 
the underneath particles motion. However, due to the 
(horizontal) roughness of the interface the effect appears 
less pronounced than in the bottom region. On the other 
hand, if the reservoir is kept open for all times, there is 
no such sudden decrease in density, and even particles at 
the interface are still quite constrained and no dense layer 
is observed in the profile shown in the previous section 

(fig- EH. 



IV. COMPACTION DYNAMICS 
A. Low density phase 

Above the jamming transition, 7 < 7 c (p R ), the ap- 
proach to equilibrium is exponentially fast, p(z,t) x 
Poc(z)+g(z)e~ t / T , as can be checked by numerically solv- 
ing Eq. 01 

When the system is in permanent contact with a reser- 
voir, the characteristic time satisfies an exact scaling 
equation [2(j 



T- 1 = 



AH 2 



(13) 



where J-(y; Pn) is a scaling function. In the special case 
of zero gravity 36], 7 = 0, particles diffuse freely from 
the reservoir until a uniform density, Poo(^) — Pri is es- 
tablished. The characteristic time of approach to equilib- 
rium can be calculated explicitly poj by linearizing Eq.El 
We find that the relaxation time for 7 = is 



AH 2 (l-p n ) 



1 



PK 
Pc 



(14) 



or equivalently JF(0; p R ) = (1 - ^/p )*/(l - p R ). Eq.Ql 
is in perfect agreement with the numerical integration of 
Eq. El 20]. As expected, the relaxation time diverges as 
Pa —> Pa- The exponent characterizing this divergence is 
<t>- 

In the presence of a gravitational field we find, by nu- 
merical integration of Eq. El that as 7 — > 7 c (p R ), the 
density of the first layer approaches p,. , p^ (0) — > p c , and 
T ~ (7c~ l)^~ 2 - Thus, the relaxation time diverges with 
exponent — 2, sec fig. El implying that the dynamics is 
faster than in the zero gravity case. Comparing with 
Eg. 1141 we see that the jamming transitions in the homo- 
geneous and inhomogeneous systems, belong to distinct 
dynamic universality classes. 



B. High density phase 

Below the jamming transition, 7 > 7 c (p R ), the density 
of the bottom layers, z < zo(p R ), is close to the critical, 
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FIG. 6: Inverse relaxation time, r _1 , for 7 = 0.1 and H = 
10 with p c = 0.88 and <f) — 3.1, corresponding to the Kob- 
Andersen model on a simple cubic lattice. The points are the 
result of numerical integration of Eq. |3 Inset shows the same 
data on the log scale. The characteristic time diverges with 
exponent <f) — 2 as the jamming transition is approached. 



p(z, t) ~ p c , and the dynamics slows down. To the lowest 
order in A(z,t) = 1 — p(z,t)/p , Eq. [21 simplifies to, 



dA(z,t) 
dt 



-7- 



dA^ 

dz 



(15) 



To solve this non-linear equation we propose a scaling 
ansatz A(z,t) — A(z/t a ). Substituting into Eq. 1151 we 
see that this form is a solution if A(z/t a ) is a power law 
with a = 1: 



A(z,t) 



l<t> t 



(16) 



Notice, that in the absence of gravity, 7 = 0, density 
relaxation is slower and characterized by a different dy- 
namical exponent, A(z,t) ~ t" 1 ^ [36|. 

Although in experiments one is usually interested in 
the bulk properties, Eq.Olcan also shed some light on how 
the upper layers z > zo compact. Using the asymptotic 
solution, Eq. 1161 in the definition of J(z,t) we find that 
for large times the particle current passing from the upper 
layers into the bulk through z — zq is 



J(z ,t) =7p 



(17) 



Since the density of the upper layers is smaller than crit- 
ical and since we are only interested in the scaling be- 
havior, it is sufficient to study the linearized version of 
Eq.El 



dp 
dt 



(Pp 

dz 2 



dp 

^ dz 



(18) 



with boundary condition: p(H,t) = p^, and J(0, t) = 
J(zo,t). The temporal Laplace transform of this linear 
equation can be easily solved yielding for the density re- 
laxation of upper layers the following expression: 



p(z, 00)- p(z,t) 



J(z ,t) 



1 



Remarkably, the time relaxation above and below z§ are 
both slow and follow a power law with different exponent. 
As expected, the dynamics in the upper layers is faster 
than in the lower ones, and its contribution to the relax- 
ation function at long times becomes negligible since <f> is 
usually larger than one. 

The asymptotic solution, Eq. 1161 is in partial agree- 
ment with the lattice-gas Monte Carlo simulation data, 
Fig. [7] The same numerical data were previously fitted 
with a four parameter logarithm law 



P{t) = Poo - 



A Pc 



l + Sln(l + t/r) ' 



(20) 



where p^ is the asymptotic packing density and B, Ap^ 
and r are adjustable parameters which also depend on x. 
The above function, first used in Ref. 0, gives a rea- 
sonable fit in the whole time window accessible to exper- 
iments; however, one can check that the long time be- 
havior is also compatible with a power law relaxation. 
Interestingly enough, something similar happens here, 
confirming that a limited time- window may not allow to 
distinguish among several regimes of slow relaxation |37| . 
One can also notice, from fig. [7| that for high values of x, 
all curves are compatible with the same exponent, while 
for small values, the exponents seems to be different. 
This may be however another artifact of the very slow 
relaxation for small x, which prevents the system from 
attaining the asymptotic regime. 

Further, we find that independently of the specific 
functional form of the fit, logarithmic or power law, the 
asymptotic packing density p^% is quite the same for 
the simulation data, Fig. |S] It turns out to be a non 
monotonic function of vibration x and displays an opti- 
mal value for the asymptotic compaction, as can be seen 
in Fig. |SJ This maximum is achieved for rather high vi- 
bration because an initial decompaction increases the free 
volume available to particles making it easier to satisfy 
the kinetic constraints and their local arrangements. It is 
clear that the specific location of the maximum depends 
on the portion chosen to measure the packing density 
(H/A in this case) however it does not affect the form of 
the plot. 

Finally, in order to stress the counterintuitive nature of 
compaction dynamics it is worth to point out its relation- 
ship with the so-called negative resistance phenomena. 
These are usually observed in a non-equilibrium station- 
ary state, where an increasing driving force leads to a 
decreasing system response (usually a particle current). 
During irreversible compaction, which is non-stationary, 
something similar happens: indeed, when vibration (the 
driving force) increases the a priori probability that a 
particle moves upward is larger; it actually turns out that 
the system compacts, i.e. the average direction of particle 
flow (the response) is preferentially downward and non- 
monotonic. This is illustrated in Fig. [5] where we plot the 
local current J(z,t) at different times t, as a function of 
the vibration strength x. One observes that at any given 
time, there is an optimum value of x at which the cur- 



7 



0.1 



0.01 

1 



0.1 



0.01 




x = 0.01 




x = 0.053 




x = 0.11 




x = 0.25 





x = 0.43 




x = 0.54 




x = 0.66 




x = 0.82 




x = 0.9 






10" 



10 2 



10" 



10° 



FIG. 7: Power law fit (normalized to unity), for large times, 
of compaction data. The fit parameters are all dependent on 
x, although a unique exponent for the large values of x is 
also consistent. The full straight line f" 1 /^ -1 ) represent the 
long-time prediction of the nonlinear diffusion equation. 
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FIG. 8: Square and circle symbols are the asymptotic pack- 
ing densities obtained within the Monte Carlo simulation with 
fixed number of particles. The two refer to the extrapolations 
using a power law fit and the logarithmic form, Eq. 1201 The 
equilibrium density is the thick solid line while the thin line 
is the bulk density evaluated from the asymptotic limit of 
Eq. I1UI and 1111 The dynamical jamming transition is lo- 
cated where all curves meet, near x c ~ 0.979. Notice that, 
from the theory, up to the region near the maximum we have 
zo > H/4 and the packing density is p c , while above it the 
point zo penetrates in the defined bulk region and the density 
deviates from p c . Inset: region near the transition. Notice 
that the simulation data agree well with the result from the 
diffusion equation, while below x c both start to deviate from 
the equilibrium curve. 



rent get its maximum. At increasing time the maximum 
moves towards higher values of x. Qualitatively similar- 
results are obtained for different z. 



C. Logarithmic compaction and Vogel-Fulcher law 

Up to now our discussion has been motivated by dy- 
namics which are characterized by a power-law vanishing 
mobility. However, there are various systems for which 
the mobility vanishes according to the Vogel-Fulcher 
law [ll|, 



r (p) = r o P exp 



ape 



P-Pc 



(21) 



Boutreux and de Gennes [38| have argued that loga- 
rithmic compaction is intimately related to this Vogel- 
Fulcher law which in turns derive from a Poisson distri- 
bution of voids in the systems. A similar conclusion can 
also be drawn from our approach. In this case the time 
evolution of the local particle density is governed by the 
equation 



dp d 



ap c 

P~ Pc 



1 dp 

1 — p dz 



IP 



(22) 
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FIG. 9: Local current J(z,t) for several times t and z = H/2, 
as a function of the vibration strength x. 



Again, proceeding as in the previous section, to the low- 
est order in A(z, t) = 1 — p(z, t)/p c , one finds that for the 
bottom layers (z < zq) 



dA(z,t) 
Ft 



-7t- exp 

oz 



A(z,t) 



(23) 
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As before we attempt a scaling ansatz. For very large 
times we find 



A(z,t) 



\n(t/z) ' 



(24) 



which is in agreement with the full numerical solution 
of Eq. [3J This relaxation is similar to the logarithmic 
compaction law, Eq. 1201 used to fit experimental data by 
Chicago group 0- 

For layers above Zq, the time evolution is still slow but 
follows a power law. Proceeding in the same way used to 
obtain Eq. [H we find 



P0,oo) - p(z,t) ~ t 1 



(25) 



Although the bulk dynamics is logarithmically slow, the 
compaction of upper layers is governed by a power law 
with exponent —1. 



V. HISTORY DEPENDENCE 

It is well known that dynamical effects in slow relaxing 
systems depend sensitively on the history of the sample 
after a quench in the high density (or low temperature) 
phase. These phenomena have been extensively studied 
by means of several experimental protocols. 

In the previous section we focused on the simplest sit- 
uation in which the sample is prepared in a random loose 
packed state and then the vibration is turned on and kept 
fixed to a given value of amplitude during the measure- 
ment. 

In order to test the system response, here we consider 
the effects of cyclic changes in the vibration amplitude, 
either continuously or suddenly, to a different value which 
corresponds to a state with a lower or higher asymptotic 
packing density. 



only appears for extremely slow driving rates. Sim- 
ilar results have been observed experimentally during 
the compaction of anisotropic granular materials like 
rods ^3 an d shearing induced compaction ■ The sys- 
tem presents a succession of irreversible branches which 
get closer and closer as the number of cycles increases 
(see Figs. ITTH and the corresponding inset). The slower 
the vibration change rate is, the smaller is the separation 
between the branches. In this way, for real systems the 
distance between the irreversible branches can become of 
the order of the measurement error, appearing as if there 
was just one reversible branch. It must be noticed that 
after cycling for a certain number of times, the density 
hardly changes from the lowest to the highest values of x, 
presenting a rather flat behavior. This can be explained 
using the asymptotic packing density, fig [H] unless the 
vibration is too high, the asymptotic bulk profile is flat. 
Even if the system enters in the high-x region, the flat- 
ness of the density will depends on the vibration change 
rate. Finally, we mention that for some choices of the 
parameters small hysteresis 39] loops also appear in the 
region of high vibration; the area of the loops is a func- 
tion of the vibration rate. We also expect that by using 
the Vogel-Fulcher mobility the cycles present a very sim- 
ilar behavior. Moreover, for high vibrations, dissipation 
effects may play a role. It would be interesting to experi- 
mentally study particles with different friction properties 
to check to what extent the cycle properties depend on 
these. In fig. ^2 we illustrate the effect of changing the 
turning point, that is, the maximum attained value of x 
before starting to decrease it. We notice that the density 
follows almost parallel paths, only that the maximum 
attained density is bigger the higher is the turning point. 



'Memory' effects 



A. 'Reversible'-irreversible cycles 

Experiments on glass beads 0, |j| have shown that un- 
der a cyclic variation of the vibration a system prepared 
in a random loose packed state first presents a branch 
during which the density increases as a function of the 
vibration (until high vibrations are attained and decom- 
paction starts). This branch is irreversible, meaning that 
when the vibration is decreased at the same rate, the sys- 
tem does not trace back the earlier evolution, but rather 
its density keeps growing as the vibration decreases. Ex- 
perimentally this second branch appears to be reversible, 
that is, the system seems to reach a stationary state in 
which any further cyclic variation of the vibration keeps 
the system always on the second branch. Along this 
branch the packing density is a decreasing function of 
the vibration amplitude (contrary to what happens in 
the irreversible compaction regime). 

Applying the same protocol to both the diffusion equa- 
tion and the lattice-gas we find that the reversible branch 



Another possible experiment devised to explore hys- 
tory dependence consists of measuring the short-time re- 
sponse of an aged system to an abrupt perturbation in 
x 0, . After evolving the system for a certain t w at 
a fixed x, the vibration amplitude is shifted by Ax until 
the time t w + At and then returned to x. As an example, 
in Fig. El we show the curves for perturbations applied 
to systems with different age, t w = 10 3 and t w = 10 4 . In 
both cases, Ax < 0, and the compaction rate increases 
while one would expect, from the long term behavior of 
the system, a slower relaxation for a smaller value of 
x. For the older system, t w = 10 4 , the first regime is 
hardly visible because the system is stiffer. These results, 
are consistent with earlier experimental pl| and theoreti- 
cal H3, |23,I13 works. For Ax > (not shown), although 
one would expect a faster compaction, for short times 
soon after the perturbation the system decompacts. The 
same behavior is also found for the diffusion equation 
when the system is in contact with the particle reser- 
voir. For large times after turning the perturbation off, 
the system resumes its normal behavior. Thus, because 
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FIG. 10: Bulk density as a function of time for a cycling 
variation of the driving rate within the nonlinear diffusion 
equation (Eq. [3J. Here the waiting time is 500 steps and 
the reservoir is again at p R = 0.1. Each branch evolves from 
5 x 10 3 steps from the minimum x value of 0.5 to the maximum 
one of 0.95. Inset: the same experiment with the gravity 
driven KA lattice-gas model. The parameters are the same 
as before, only that the waiting time is 5000 MCS and each 
branch has a total duration of 2 x 10 4 MCS. In both cases, 
the distance between the branches decreases as the number 
of cycles increases. 



FIG. 12: Short-term memory experiment performed with the 
gravity driven KA model. The thick line shows the unper- 
turbed evolution. The system is perturbed at t w = 10 3 and 
10 4 after evolving with x = 0.6 (solid normal lines). In the 
interval t w < t < t w + At, with At — 10 4 , the vibration 
changes to 0.3. Notice that in the case where the system is 
perturbed earlier, the system first increases its density (the 
short-term memory effect as described in the text), resuming 
its expected behavior after some time. In the case of a later 
perturbation, this first regime is not noticeable. 
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FIG. 11: The same as the inset of fig. 1101 only that x is re- 
versed at different maximum values. Notice that in this case, 
different densities are attained. 



of the transient nature of the response, it appears as a 
short-term memory. Moreover, after the initial anoma- 
lous behavior, the perturbed and non perturbed curve 
crosses and the system start to evolve with the rate ex- 
pected from long term behavior. 

Besides short-term memory effects, long-term memory 
is also present. To see this, we repeat the above exper- 
iments with the difference that the perturbation is ap- 
plied at a much larger t w (10 5 in this case). Moreover, 
the perturbation is kept on also for a longer period. After 
being turned off, the system vibration is returned to its 



previous values. As can be seen in figure ITSI the system 
evolves much less during the time interval At when it is 
perturbed: the h cm roughly follows a plateau. This is 
expected since the relaxation rate decreases for smaller 
values of x. Moreover, when turning the perturbation 
off, the system returns to a point very close to the one 
where it was before being perturbed. This is more clear 
in the inset of the same figure, where the perturbed data 
for t > t w + At are shifted by At and seems, as a first ap- 
proximation, to collapse on the unperturbed curve, show- 
ing that the system keeps, to some extent, memory of its 
state before the perturbation even if being perturbed for 
a long time. 

As is common is systems in a non-stationary state, the 
response of the system to a perturbation in the vibration 
clearly depends on its age. The larger is t w , the greater is 
the compaction achieved by the bulk and less responsive 
the system becomes. For small t w , the bulk is still very 
sensitive to the perturbations, the amount of free volume 
is considerable and as soon as the vibration is lowered, 
the particles get closer and the density increases very 
fast, even faster than one would expect from the knowl- 
edge of the asymptotic behavior. By increasing, instead 
of decreasing, the vibration, the opposite behavior is ob- 
served. Thus, short-term memory is related to pertur- 
bations at early times. However, as t w increases, the 
amount of free volume in the bulk decreases and most of 
the instantaneous response, when perturbed, comes from 
the interface, that has a fast dynamics. After this strong 
response due to the interface (that is only seem in global 
measures like the height of the center of mass), the system 
continues at a much smaller pace, corresponding to the 
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A. Mean-square displacement 

In our approach the aging effects are best studied by 
considering the two-times mean square displacement of 
particles, B(t,t w ), which for the 3D lattice-gas is defined 
as 



3 N 

s (*< *») = E E + o - r S(*-)] 2 ) 



(26) 



o=l fc=l 



where rjjj(i) are the coordinates (a = 1, 2, 3) of the parti- 
cles k at times t. In the continuous ID diffusion model, if 
t is sufficiently larger than t w , the mean-square displace- 
ment at height z can be written as 



FIG. 13: Long-term memory experiment. The evolution of 
the height of the center of mass (or, equivalently, the potential 
energy) is plotted as a function of time. Notice that if the 
system has evolved enough time (t w — 10 s in this case) the 
further evolution after the perturbation is quite small and the 
system returns to its previous state when the perturbation 
is turned off. In the inset, the perturbed data for t > t w 
are shifted, t — » t — 10 5 showing that they collapses on the 
unperturbed data. Averages are over 300 samples. 



expected evolution at the new vibration value. Thus, in 
the interval At while the perturbation is on, the amount 
of change is due to the bulk evolution and is smaller the 
larger is t w (for example, for the waiting times that we 
probed, from 10 3 until 10 5 , the bulk always evolved). 
As a conclusion, "long-term memory" effects can only be 
seen as an approximation, with a very simple explanation 
as also pointed out in |24| . 

As remarked in Refs. jljj and ^lj ( m the context of 
glasses) , a rate equation for a single macroscopic variable 
(e.g., the free volume) would not be able to account for 
the complexities of the memory effects. Nevertheless, it 
must be emphasized that the nonlinear diffusion equa- 
tion studied here is explicitly dependent on the spatial 
dimension, and is this heterogeneous profile that encodes 
the additional information responsible for the effects dis- 
cussed here. However, more complex memory effects are 
observed in systems without gravity like glasses and spin- 
glasses, and is unlikely that this equation would be able 
to account for them, since probably they involve effects 
due to the interaction. In order to capture the complexity 
of these effects, this equation would have to be properly 
generalized. 



VI. PHYSICAL AGING UNDER GRAVITY 

We now turn to the discussion of physical aging phe- 
nomena under gravity as they appear in the two-time 
scaling behavior after a sudden quench into the high 
packing density phase. 



B z (t,t w ) 



ds r [p(z, s)] 



(27) 



It is worth to recall that in zero gravity one finds a simple 
aging 36] . For a power-law diffusion and for the lower 
layers z < zo(Pr) we find, to leading order in t and t w , a 
two-time scaling of the form 



B z (t, t w ) ~ t 



- t 



(28) 



with an exponent /i = 4>/{4> — 1). Since, usually <f> > 1, 
this corresponds to a super-aging regime, fi > 1. This 
means that the effective structural relaxation time grows 
as t!£, what is faster than the age of the system, t w . It 
also means that there is a microscopic time scale that 
starts to become relevant, differently from what happens 
with simple aging. For the upper layers, z > z (p R ) one 
obtains the same two-time scaling behavior but this time 
the exponent is [i = 4> 2 /(</> — 1). However, for closed 
systems at not so large vibrations, the contribution from 
the layers above zo is small since most of the particles 
are at z < zq. 

A similar super-aging behavior has been observed in 
the simulation of the gravity-driven KA model ^3 where 
the agreement is rather good if the vibration is not too 
low. For example, we show in fig. 1141 the mean square 
displacement for the case x = 0.4 along with fig. El 
where these curves were collapsed onto a master curve 
following the above scaling. The super-aging exponent 
obtained from the data collapse is fj, = 1.48 which is 
in very good agreement with the theoretical prediction 
(j, = <f>/(4> — 1) — 1.476. For smaller vibrations, the ex- 
ponents become smaller since the time accessible to the 
measurement is probably not enough to reach the asymp- 
totic regime (where the approximation is valid) , thus dif- 
fering from the theoretical value. 

For a Vogel-Fulcher diffusion instead we find, to the 
leading order in t and t w , and for z < zo(p R ): 



B,(t,tw)~log( ^ 



(29) 



that is a simple aging scenario, whereas for the upper 
layers, z > zo(p R ), one obtains (to the leading order in 
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FIG. 14: Mean squared displacement B(t,t v ) for waiting 
times t w = 2 10+k (k = 1,2,...) and vibration amplitude x = 
0.4. 
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FIG. 15: Mean squared displacement B(t,t w ) for several 
values of t w and x = 0.4 as a function of t w ~'* — t The 
data collapse is obtained for n = 1.48. 

t and £ w ) that the two-time scaling does not depend on 
t/t w but rather 

B,(t,«w) ~e- Q ^*" -e~ ap ^. (30) 

B. Triangle relation 

It is interesting to observe the different behavior of 
Eas.1291 and 1281 at finite waiting times t w . In the case of 
simple aging, ]im.t-nx B(t,t w ) = oo, i.e., a weak ergod- 
icity scenario [43j; while for super-aging, a finite limit is 
obtained (which, however, vanishes as f w — * oo). The 
manner in which time-translation invariance is violated 
is, however, similar. Indeed, if we consider the "triangle 
relation", B(*i,t 3 ) = f[B(ti,t 2 ), B(t 2 ,t 3 )], where the 
times t\, t 2 , and t% are in increasing order, it is straight- 
forward to check that f(x,y) = x + y in both Vogel- 
Fulcher and power law cases, implying that displacements 



over non-overlapping time intervals are statistically inde- 
pendent. This feature does not hold in the presence of 
activated aging for which B(t,t w ) ~ logi/logi w |43|. In 
the gravity driven KA model the triangle relation is not 
obeyed at short times but becomes asymptotically valid 
at longer times as it can be seen in fig. ^3 5jJ . 
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h 

FIG. 16: Check of the triangle relation for the mean squared 
displacement B(t,t') for x — 0.1. Different symbols stand for 
different pairs {t2,tz) as indicated in the figure. 



VII. DYNAMIC HETEROGENEITIES 

In purely kinetic models, the absence of an increas- 
ing static correlation lenght on approaching the dynamic 
arrest poses the question of whether the increasing relax- 
ation times can be related to a diverging dynamic cor- 
relation lenght. This has been associated to the pres- 
ence of dynamical heterogeneities in glasses (for a review 
see [25|, |26|, |33 and references therein). If the glass/dense 
granular analogy holds true then one would expect that 
the role played by dynamical heterogeneities in slow gran- 
ular compaction should be similar to that observed in 
glassy dynamics. However, the role played by these 
structures and the associated lenghts, on the dynamics 
of granular and colloidal systems, is yet to be under- 
stood. Several measures for quantifying the spatial het- 
erogeneities have been introduced for kinetic models . 
In particular, this issue was recently investigated in the 
KA model |44j | without gravity using a fourth-order cor- 
relation function. Here we study the role played by these 
structures in the non-zero gravity case. 

In figure H7I we plot the dynamical nonlinear response 

Xa(z, t) = N ({q 2 (z, t)> - (q{z, i)) 2 ) (31) 

where N is the number of involved sites in the computa- 
tion, q(z,t) = C(z,t)/C(z,0) and 

C{z ^ *) = 4 H ^(*H(°) " °) ( 32 ) 
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where i runs over all sites in the z, z — 1 and z + 1 layers. 
Consistently with the theoretical expectation [3^, the 
long time limit of \i converges to unit, Xi( z i 00 ) = 1 [H3- 
We also verified that this asymptotic behavior is valid 
in absence of gravity at variance with the results of 
Ref. 0| W^i- Analogously to what happens in the KA 
model without gravity and in other glassy systems, the 
peak is shifted to higher times and gets larger as the den- 
sity increases (the lower is z, the greater is the density). 
In the inset of fig. El we show that both the position 
and height of the peak grow as power laws as the density 
of the corresponding height approaches p c . Interestingly, 
Xa only depends on the local density: for example, in 
fig. 1171 is shown that two curves corresponding to differ- 
ent z and x, but having almost the same density, are 
coincident within the numerical accuracy. 
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FIG. 17: Dynamical response, eq. 1311 as a function of time 
(in MCS) for different vibrations: x=0.92 (filled symbols) and 
0.94 (empty symbols). Different symbols stand for different 
heights: z = 5, 10 and 15 (square, circle, triangle, respec- 
tively). The line is the asymptotic X4 = 1 behavior. Notice 
in the figure the presence of two very close curves: they cor- 
respond to different vibrations and heights, but their density 
is the same within the numerical accuracy. Inset: location 
(circles) and height (triangles) of the maximum of \i- Both 
diverge as power laws with approximate exponents 3 and 0.6, 
respectively. 



VIII. CONCLUSIONS 

We have investigated some aspects of slow gran ular 
dynamics inspired by kinetic lattice-gas models |15l l31j . 
The key ingredient of these models is a free-volume re- 
striction implemented by a purely kinetic constraint. No 
interaction between the particles is assumed beyond the 
hard core exclusion. The thermodynamics of the model 
is completely trivial and all its interesting features are 
purely dynamical. A macroscopic transport equation was 
written and studied, allowing us to predict the specific 
location of a jamming transition and to analyze the be- 
havior of the system in its vicinity. 



As noted in the earlier work |2fJ and detailed here, 
the time evolution of particle density in a gravity driven 
lattice gas is completely controlled by the mobility of 
the corresponding gravitationless system. It was shown 
that for a power law mobility, the bulk density relaxes 
as a power law. On the other hand, when the parti- 
cle diffusion decreases exponentially accordingly to the 
Vogel-Fulcher law, a logarithmically slow compaction is 
found. Due to the finite time window available, the data 
on granular systems from the Monte Carlo simulation, 
as well as from experiments, is consistent with both the 
power law and logarithmic relaxation. Furthermore, the 
issue of whether the asymptotic density is a monotonic 
function of the vibration amplitude, can only be solved 
by performing simulations over a much larger time win- 
dow |37| . For the maximum time achieved here, the 
asymptotic packing density seems to depend on the vi- 
bration strength in a non trivial way. 

The similarities between granular and glassy systems 
have been stressed many times in recent years. Here we 
extended this discussion by showing that the behavior 
of dynamical heterogeneities is quite similar to the ones 
present in systems without gravity. We showed that for 
different vibrations the global non homogeneity induced 
by gravity (that is, the density profile, different for each 
vibration) does not affect the local character of these 
quantities that only depends on the local density: the 
heterogeneities present in a horizontal layer depends only 
on this layer density. This result could also be relevant 
to the investigation of slow sedimentation of colloids. 

Other signatures of glassy behavior are also present, 
like aging, reversible-irreversible cycles and memory ef- 
fects. In particular, short and long-term memory effects 
are simpler than their glassy counterpart and have al- 
ready been described in terms of the density profile prop- 
erties. Irreversible branches are obtained when cycling 
the vibration amplitude, approaching a quasi-reversible 
branch when the rate is slow enough or after cycling 
many times. We also remark the analogy between short- 
term memory effects and the cycling experiments. In the 
decreasing- X part of the cycles, the density increases as 
for smaller x. This is at variance with what one would 
expect from the long time behavior of the asymptotic 
packing density, exactly in the same way as the system 
behaves in short-term memory experiments. The differ- 
ence is that in one case the change in vibration is discon- 
tinuous while in the other, it is done at a small rate. 

The analytical approach presented here has some lim- 
itations. Among the features seen in the simulation, 
which are not properly described by the theory, are: exis- 
tence of a dense layer between the bulk and the interface 
and the initial state dependent oscillations at the bot- 
tom of the sample. These may be a direct consequence 
of the local density approximation used for the mobil- 
ity. A weighted density approximation might be able to 
account for some of these features. Moreover, since the 
equation is deterministic, fluctuation dependent quanti- 
ties (e.g., Xi(t)) are not captured by the formalism and 
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noise has to be included |45j . 

Another interesting question concerns the role of fric- 
tion. In realistic systems, the energy injected by stirring, 
shearing, tapping etc. is dissipated by the inelastic colli- 
sions among the particles, creating a gradient of granular 
temperature 0, E3> ■ Here we assume that the sys- 
tem is in contact with a thermal bath, a situation which 
is quite different. In highly packed granular systems how- 
ever the steric hindrance is far more important than the 
energy dissipation, as one is able to reproduce several 
features of compaction dynamics by ignoring the specific 
mechanism of energy injection/dissipation. 

Finally, we mention that slow relaxation in the KA 
model has also been analyzed in terms of the dynamically 
available volume, expressed by holes density v [4j|. A 
hole is defined as an empty site where a neighbouring par- 
ticle satisfies the kinetic constraints and may jump in. It 
has been recently suggested that dynamical arrest in dif- 
ferent systems, such as glasses, colloids, gels, compressed 
emulsion and granulars, has a universal character when 
described in terms of v. In particular, as one approaches 
the transition, the diffusion coefficient goes to zero as a 
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